Climate change multi-model projections in CMIP6 scenarios in Central Hokkaido, Japan

Simulation of future climate changes, especially temperature and rainfall, is critical for water resource management, disaster mitigation, and agricultural development. Based on the category-wise indicator method, two preferred Global Climate Models (GCMs) for the Ishikari River basin (IRB), the socio-economic center of Hokkaido, Japan, were examined from the newly released Coupled Model Intercomparison Project Phase 6 (CMIP6). Climatic variables (maximum/minimum temperature and precipitation) were projected by the Statistical DownScaling Model (SDSM) under all shared socioeconomic pathway-representative concentration pathway (SSP-RCP) scenarios (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP4-3.4, SSP4-6.0, SSP5-3.4OS, and SSP5-8.5) in two phases: 2040–2069 (2040s) and 2070–2099 (2070s), with the period of 1985–2014 as the baseline. Predictors of SDSM were derived from CMIP6 GCMs and the reanalysis dataset NOAA-CIRES-DOE 20th Century Reanalysis V3 (20CRv3). Results showed that CMIP6 GCMs had a significant correlation with temperature measurements, but could not represent precipitation features in the IRB. The constructed SDSM could capture the characteristics of temperature and precipitation during the calibration (1985–1999) and validation (2000–2014) phases, respectively. The selected GCMs (MIROC6 and MRI-ESM-2.0) generated higher temperature and less rainfall in the forthcoming phases. The SSP-RCP scenarios had an apparent influence on temperature and precipitation. High-emission scenarios (i.e., SSP5-8.5) would project a higher temperature and lower rainfall than the low-emission scenarios (e.g., SSP1-1.9). Spatial–temporal analysis indicated that the northern part of the IRB is more likely to become warmer with heavier precipitation than the southern part in the future. Higher temperature and lower rainfall were projected throughout the late twenty-first century (2070s) than the mid-century (2040s) in the IRB. The findings of this study could be further used to predict the hydrological cycle and assess the ecosystem's sustainability.


Materials and methods
Study area. The Ishikari River basin (IRB, 42°41′9.6″N-44°47′8.9″N, 140°59′33″E-143°10′47.8″E) was situated in the Mid-western Hokkaido, Japan, with an area of 14,330 km 2 (Fig. 1). The Ishikari River originates from Mt. Ishikari and flows westward to the Sea of Japan. Toyohira, Tobetsu, Chitose, and Yubari are its major tributaries 34 . 52% of the population of Hokkaido live in the IRB. The IRB is an important economic, agricultural, industrial, and cultural center of Hokkaido, and is also the seat of Sapporo and Asahikawa, the largest two cities of Hokkaido. The IRB is dominantly controlled by the hot-summer subtype and hemiboreal climate. According to the climate records (1985-2014) of 13 meteorological stations across the IRB provided by the Japan Meteorological Agency (JMA, http:// www. jma. go. jp), the annual average maximum and minimum temperatures are 10.27-12.74 °C and 0.34-5.67 °C, respectively. The annual precipitation is about 1007-1610 mm. The rainy season is generally from August to September. The snowfall period is from mid-December to late March of the following year, with an average annual maximum snow depth of 35 cm. Hydrologic peaks occur during the snow-melt period (March to May). Data collection. Three types of meteorological datasets were employed in this study, including observed historical data, reanalysis data, and GCM data. IRB covered 13 meteorological stations (as shown in Fig. 1). The historical meteorological data across 30 years period of 1985-2014 was composed of daily maximum air temperature at 2 m (tasmax), daily minimum air temperature at 2 m (tasmin), and daily precipitation (pr) of each meteorological station, which can be accessed in the JMA.
In this study, the reanalysis dataset selected the most recent version of reanalysis from the Twentieth Century Reanalysis (20CR) Project, which is funded by the National Oceanic and Atmospheric Administration (NOAA), the Cooperative Institute for Research in Environmental Sciences (CIRES), and the U.S. Department of Energy (DOE), NOAA-CIRES-DOE 20th Century Reanalysis V3 (20CRv3). The 20CRv3 contained objectively-analyzed 4-dimensional weather maps and their uncertainties 35 . The 20CRv3 covered the spatial resolution of 1.0-degree latitude × 1.0-degree longitude global grid (360 × 181). Daily atmospheric variables (such as humidity, precipitation flux, and geopotential height) of 20CRv3, spanning 1985 to 2014, were readily accessible at the NOAA Physical Sciences Laboratory (PSL, https:// www. psl. noaa. gov/).
Monthly and daily GCM datasets under the first variant level "r1i1p1f1" were obtained from the Coupled Model Intercomparison Project Phase 6 (CMIP6, https:// esgf-node. llnl. gov/ search/ cmip6/). According to the www.nature.com/scientificreports/ availability of predictors, 17 GCMs were selected in this study, as shown in Table 1. All these predictors were interpolated onto a 1° × 1° grid to match the 20CRv3. The GCMs-derived monthly variables (tasmax, tasmin and pr) in the period of 1985-2014 were applied for the selection of GCMs. The daily atmospheric predictors (corresponding to 20CRv3), spanning 1985-2100, were applied in statistical downscaling analysis to generate and project future climatic variables under eight Shared Socioeconomic Pathways-Representative Concentration Pathways (SSP-RCPs) scenarios (SSP1-1.9, SSP1-2.6, SSP4-3.4, Figure 1. Locations of the Ishikari River basin (IRB) and its meteorological stations. Note: This figure was generated by ArcMap 10.2 (https:// suppo rt. esri. com/ en/ produ cts/ deskt op/ arcgis-deskt op/ arcmap/ 10-2-2). Shapefiles are available from the Japanese Geographical Survey Institute (JGSI, http:// nlftp. mlit. go. jp). www.nature.com/scientificreports/ SSP2-4.5, SSP4-6.0, SSP3-7.0, SSP5-3.4OS, and SSP5-8.5). CMIP6 used a matrix framework that combined two determinants of emission scenarios (e.g., RCPs) and a diverse range of socioeconomic assumptions, namely the so-called Shared Socioeconomic Pathways (SSPs) scenarios, to force climate models 17,36,37 , which makes future scenarios more reasonable. The growth of civilization and natural systems at the national and regional levels in the twenty-first century provides a foundation for the formation of SSPs 38 . Five narratives of SSP (SSP1, SSP2, SSP3, SSP4, and SSP5) scenarios were employed to picture the potential challenges brought about by the variations in global and regional evolution across time 39,40 . Numbers are consistent with low-to-high concerns in mitigation and adaptation in the futuristic society, named SSP1 ("taking the green road"), SSP2 ("a middle of the road"), SSP3 ("regional rivalry-a rocky road"), SSP4 ("inequality-a road divided"), and SSP5 ("fossil-fueled development"), respectively 41-43 . Evaluation of GCMs performance. GCMs are extensively applied to simulate past climates and produce future climatic variables 44 . However, there is a significant uncertainty in estimating regional applications of GCMs due to the difference in each GCM, such as resolution (fine or coarse), climatic response mechanism (aerosols, circulations of land, ocean, and atmosphere), and spatial-temporal scales 45 . Hence, there is an urgent need to analyze each chosen GCM to minimize the uncertainties when applying them in specific areas. Evaluating the performance of GCMs simulation is generally to compare them with reanalysis or observed climatic data. Lots of indicators have been employed by various researchers in climate modelling. Raju and Kumar reviewed more than hundreds of works on climate models to study which are the best GCMs 46 . They recommended using category-wise indicators when evaluating GCMs, such as error and correlation coefficient. The Taylor diagram can compare simulations (model) with measurements using the correlation coefficient, root-mean-square difference, and standard deviations to graphically assess these qualities 47 . The high correlation and few errors represent that selected GCMs are suitable for the local climate system. Therefore, the Taylor diagram was used in this work to compare CMIP6 outputs with regional data and further to identify the best CMIP6 GCMs for modeling temperature and precipitation across the IRB.  23 . SDSM can set up statistical relationships between large-scale predictors and regional-scale climatic conditions (e.g., temperature and precipitation) using a combination of multiple linear regressions. If these correlations remain true as in prospective, they may be utilized to acquire downconverted regional features in certain coming phases by forcing the interactions with GCM-derived predictors through the stochastic weather generator. There are two different progressions in each sub-model, unconditional and conditional. Temperature does not need to be transformed and directly generated in the unconditional pattern, which exhibits a linear relationship between the predictors and predictand (e.g., individual wind speeds can be used to calculate regional airflow parameters). Precipitation should be reformed by the fourth root and then simulated in the conditional pattern, which is an intermediate process between regional forcing and local climatic conditions. For instance, local precipitation is determined by the occurrence of rainy days, while the latter is determined by regional-scale predictors (such as moisture and atmospheric pressure). The wet criterion for daily rainfall was chosen at 1.0 mm in this study, as it is commonly employed for statistical downscaling.
Model process. SDSM performs five key steps, from variables selection, calibrating and validating model, to weather generation and future climate scenarios projection 48 . Screening variables between predictand (such as maximum temperature, minimum temperature, evaporation, as well as precipitation on a local scale) and predictor (large-scale atmospheric conditions) is a core of the statistical downscaling process. SDSM combines the correlation matrix, partial correlation, P value, histograms, and scatter plots, which can help users to find the best predictors. Peng et al. provided extensive explanations of each SDSM procedure 49 .
In this study, predictors were constructed by the reanalysis datasets (20CRv3) and CMIP6 GCMs to reproduce ensembles of present climate data in SDSM. Commonly used predictors are normalized and obtained as predictor datasets. In this study, predictors are comprised of mean temperature at 2 m (temp), mean sea level pressure (mslp), total precipitation (prcp), surface downwelling longwave flux in air (rlds) and surface downwelling shortwave flux in air (rsds) in monolevel, specific humidity (#_shum), relative humidity (#_rhum), geopotential height (_p), geostrophic air flow velocity (#_f), vorticity (#_z), zonal velocity component (#_u), meridional velocity component (#_v), divergence (#_zh), and wind direction (#_th) under the pressure level. Table 2 lists the most suitable predictors for observed predictands in 13 meteorological stations across the IRB. Simulation of precipitation needs more predictors other than temperature.
SDSM evaluation. Generally speaking, the model's performance is mainly based on the selection step for predictors to predictands in SDSM. Even though there are statistical or graphical ways, such as correlation matrix and P value during the screening process, to identify the most accurate predictors, applying statistical parameters in the evaluation process is necessary to avoid uncertainties in SDSM. The ability of modeling outputs from SDSM was obtained by determining the Nash and Sutcliffe efficiency [NSE, Eq. The magnitude of NSE is computed using the formula below: www.nature.com/scientificreports/ The magnitude of R 2 is calculated using the following equation: The magnitude of RMSE is computed using the formula below: The magnitude of Pbias is calculated using the following equation: where X oi is the observed predictand on day i, X mi is the modeling outcome on day i, X oi is the average measured value during the study period, and n is the total number of the observed data. The NSE illustrates how well the observed and simulated data suit the 1:1 line. Both R 2 and RMSE are indices of quality of fit, whereas Pbias reveals the model's tendency to over-or under-estimated with respect to the observed data. The model performs well when R 2 and NSE values are close to one, and the lower the RMSE and absolute value of Pbias are, the tighter the modeled and measured magnitudes are 50 .

Results
Selection of GCMs. When assessing GCMs, it is vital to compare GCM outputs with observed records; otherwise, even previous robust predictions may not provide skilled future projections 51 Table 2. List of selected predictors in the Ishikari River basin. Mean temperature at 2 m (temp), mean sea level pressure (mslp), total precipitation (prcp), surface downwelling longwave flux in air (rlds) and surface downwelling shortwave flux in air (rsds) in monolevel, specific humidity (#_shum), relative humidity (#_ rhum), geopotential height (_p), geostrophic air flow velocity (#_f), vorticity (#_z), zonal velocity component   SDSM downscaling. In this study, the span of 1985-1999 was considered as the calibrated phase, and the validated phase was from 2000 to 2014. Table 3 lists the values of four evaluation metrics (such as R 2 , RMSE, NSE, and Pbais) at 13 meteorological stations when calibrating and validating maximum/minimum air temperature and precipitation on a monthly scale. The R 2 and NSE values of each meteorological station's maximum or minimum air temperatures throughout the simulation phase were nearly equal to 1.000. In terms of maximum air temperature, the values of RMSE in these two modeling phases ranged within 0.035-1.065 and 0.089-0.964, respectively. The Pbias values were less than 0.400% in the calibration phase, with the exception at the Horokanai  www.nature.com/scientificreports/ the calibration phase performed better than that in the validation phase. Meanwhile, the outcomes demonstrated that SDSM outperforms precipitation in modeling tasmax and tasmin. Precipitation is always difficult to simulate due to its high dynamic properties 52 . Generally, choosing an effective mix of predictors for SDSM is challenging owing to influencing variables such as dry and wet period durations, local microclimates, and terrains, which could not be fully covered in reanalysis datasets 53 . In this study, results in downscaling temperature and precipitation demonstrated that predictors from the 20CRv3 dataset are able to reflect attributes on a local scale. In total, SDSM performs effectively in simulating the climatic variables during both calibration and validation phases in the IRB.  Annual variations in precipitation showed more discrepancies not only at different meteorological stations but also under different SSP-RCPs scenarios from two GCMs (Fig. 7). For example, at the Pippu, Kamikawa, Takikawa, and Ashibetsu stations, SSP-RCPs scenarios from MIROC6 generated higher annual precipitation than those in MRI-ESM-2.0. But at other stations such as Horokanai, Scenario SSP1-1.9 from MRI-ESM-2.0 could predict more precipitation than that in MRIOC6. As shown in Fig. 7, annual precipitation was expected www.nature.com/scientificreports/   (Fig. 9). The emission scenario had a greater impact on temperature and precipitation projection than the socioeconomic scenario. The increasing range is from the high-emission (SSP5-8.5) to the low-emission case (SSP1-1.9). Those results correspond to the original setting of each scenario of CMIP6. In direct contrast, precipitation exhibited opposite variations under SSP-RCP scenarios compared to tasmax and tasmin. Under MRIOC6, SSP4-3.4 (24.1 mm) produced the greatest increase of rainfall and SSP3-7.0 (− 38.7 mm) made the least in 2040s, SSP1-1.9 (30.3 mm) generated the highest amount of rainfall and SSP5-8.5 (− 136.5 mm) made the greatest decrease in 2070s. Under MIR-ESM-2.0, SSP1-1.9 (− 15.0 mm) generated the greatest amount of rainfall and SSP5-8.5 (− 124.8 mm) made the least in 2040s, SSP1-2.6 (− 35.0 mm) produced the greatest amount of rainfall and SSP4-6.0 (− 128.5 mm) made the least in 2070s (Fig. 8). In the comparison of each SSP-RCPs scenario, SSP1-1.9 always produced the highest amount of rainfall, with the percent changes of − 0.3% (mean value of two GCMs) in the 2040s and − 1.1% in the 2070s, and SSP5-8.5 made the least, with the percent changes of − 6.3% (2040s) and − 10.9% (2070s) (Fig. 9). In addition, MIROC6 could support a wider range of variations in maximum and minimum temperatures than MRI-ESM-2.0. The percent changes under MIROC6 were 16.9-51.5% for tasmax and 77.5-220.3% for tasmin, respectively. In terms of MRI-ESM-2.0, variations were from 18.2 to 46.8% for tasmax and from 77.  www.nature.com/scientificreports/ The distribution of average changes in maximum temperature, minimum temperature, as well as precipitation of the IRB, during the 2040s and 2070s under the MIROC6 and MRI-ESM-2.0, are displayed in Figs. 10, 11, and 12, respectively. MIROC6 and MRI-ESM-2.0 both show the higher air temperature and less precipitation (Figs. 10, 11, and 12). Spatially, when it came to the distribution of tasmax and tasmin, both products showed a similar trend, with a bigger rise in the northern part and a smaller rise in the southern part of the IRB. According to meteorological data, the northern part was colder than the southern part of the study area, but CMIP6 GCMs anticipated that the northern part may exhibit a stronger warming trend in the future. The precipitation across the IRB showed a decreasing trend under MRI-ESM-2.0. However, a large increase in precipitation produced by MRIOC6 was found in the mid-northern part of the IRB, which was displayed as an absolute difference from MRI-ESM-2.0. Temporally, higher temperature and less precipitation were projected during the late twentyfirst century (2070s) than the mid-century (2040s) in the IRB. The distribution of climatic variable changes was strongly affected by emissions. For all SSP-RCP scenarios assessment, the distribution of climatic variables change was strongly affected by emissions (Supplementary Figs. 1-3). Higher emissions are associated with higher temperatures and less precipitation under all of the climatic scenarios in this study.

Discussion
In this study, we compared climatic variables (air maximum/minimum temperature and precipitation) from 17 available CMIP6 GCMs at 13 meteorological stations in the IRB. The performances of 17 GCMs were evaluated by Taylor diagrams with two evaluation metrics (correlation coefficient and root-mean-square difference). Most CMIP6 GCMs were able to access temperature measurements with high correlation coefficients. But there is not a single GCM could well reproduce the precipitation characteristics across the IRB. Our results showed only two preferred GCMs, i.e., MIROC6 and MRI-ESM-2.0 showed the best adaptability in temperature and precipitation across the target region.
Accordingly, in order to generate daily maximum/minimum temperature and daily precipitation in two future phases, the 2040s (2040-2069) and 2070s (2070-2099), a statistical downscaling model was established based on the 20CRv3 reanalysis datasets and CMIP6 GCMs-derived predictors, with respect to observed climate during   14 . That is because GCMs are incapable of modeling precipitation with high accuracy, e.g., on the local station scale, and the GCM outputs of precipitation are affected by topographical factors and regional climatic forcing 23 . On the other hand, the downscaling process of precipitation further propagates this error 57 . It should be more cautious when downscaling precipitation with SDSM 58 . To sum up, FGOALS-g3, CanESM5, MIROC6, and MRI-ESM2-0 have the best predictive ability of maximum temperature, MRI-ESM2-0 has the best predictive ability of minimum temperature, and MIROC6 has the best predictive ability of precipitation. Given the availability of each GCM with adequate predictors for SDSM simulating, two GCMs, i.e., MIROC6 and MRI-ESM2-0, from CMIP6 were chosen to generate future climatic variables under various SSP-RCPs scenarios.  (Fig. 8). The average temperature of two GCMs was anticipated to rise by 1.89-3.10 °C (2040s) and 1.91-5.02 °C (2070s) under all scenarios. The latest generation simulations, CMIP6, continue to paint a picture of more frequent and intense high temperatures on a wide scale 59 . It's hardly surprising that projected climate enhanced warming during 2040-2099 under all SSP-RCP scenarios. On a global level, the annual mean temperature under CMIP6 shows a more sensitive warming trend at higher latitudes 60 , which indicates that the warming effect in the north region of IRB is greater than that in the southern (Figs. 10 and 11).
Future annual rainfall, generated from MIROC6 and MRI-ESM-2.0, broadly decreased in the projected periods (2040s: − 1.0% and − 6.3%; 2070s: − 3.4% and − 7.0%, respectively) (Fig. 9). When the air temperature rises, the increasing evaporation may be able to remove more water, resulting in more frequent and intense storms or less precipitation and increased risk of drought 61 . In the example of the Ishikari River basin, it may explain that the climate with a greater warming effect from downscaled GCMs derived may cause the less total amount of precipitation over land. Similar findings have appeared in the Duan et al. which revealed the decreased precipitation amount under CMIP3 GCMs in the upper Ishikari River basin 62 , this is consistent with our findings. Meanwhile, the distribution of rainfall under CMIP6 GCMs shows discrepancies in variable rate between the north and south regions in Ishikari River basin. A number of studies indicate that warming will increase www.nature.com/scientificreports/ precipitation in equatorial and high-latitude regions 63,64 . Generally, our finding suggested that more amounts of precipitation would appear in the northern than southern of the IRB (Fig. 12). However, uncertainties in climate models are inevitably subjected to prove the limiting factors, contributing to model structures 65,66 , downscaling approaches 67,68 , and local forcing, even for very broad regional averages 69 . When further employing the GCMs-derived datasets, how to reduce this uncertainty is not unavoidably overlooked. Different climate models are constructed based on sophisticated amalgamations of various patterns and assumptions, which would lead to large uncertainties in this work of climate projections 70 . But several studies have exemplified that the multi-model ensemble could be an efficient method to lower model-specific uncertainty 71 . Here 17 available CMIP6 GCMs and eight scenarios were employed to construct the future climate (temperature and precipitation) in the IRB. CMIP6 GCMs were selected by evaluating model applicability with several metrics. Preferred multi-GCMs and SSP-RCP scenarios could provide a wider variation range of climatic variables to eliminate further uncertainties when projecting future climate change scenarios 72 . The findings of our work could support the future projections of essential climatic variables under all the SSP-RCP scenarios. The uncertainty envelope was also largely shaped by the downscaling techniques. There should be caution if only one approach is used to reduce uncertainty in downscaling progress. Particularly, using the regression-based methods (such as the SDSM application in this study) would raise the concern of the anomalous downscaled future temperatures increasing 73 . Besides forcings of model and downscaling, this, specific geographic location of the IRB, in particular, makes precise precipitation distribution and amount predictions difficult 74 . Due to the topography of the coast and the complexity of orographic effects (Mashike Mountains and Taisetsu Mountains). During the Asian winter monsoon, sea-effect snowfall generated over the Sea of Japan affects the Ishikari River basin 75,76 . In addition, several efforts on climate modeling have pointed out that global warming causes more intense precipitation but reduced precipitation frequency 63 . Since precipitation is such complex combinations of different physical processes, the volume of precipitation is determined not only by the severity but also by the frequency 77 . In this study, we just estimated the amount of precipitation, lacking other characteristics of precipitation (e.g., intensity and frequency), further study is needed to properly elucidate the influence of global warming on the internal variability of the geophysical system, and then to fully understand how climate change drives the mechanisms that result in extreme events 71 . Since the climate with warming and less rainfall would pose a significant challenge to ecosystems of IRB, some studies show rising temperature can increase crop yield 78 . Groundwater flow is also affected by precipitation when less streamflow caused by rainfall decreasing 79 . Those fundamental datasets would be dedicated to further climate change related work, such as modeling hydrological processes, simulating crops www.nature.com/scientificreports/ growth and yield, evaluating ecosystem services, and so on. Those potential situations should be considered in further studies when developing models.

Conclusion
In this study, reconstructed SDSM was successfully developed using the datasets from the 20CRv3 and two preferred GCMs (MIROC6 and MRI-ESM-2.0). Then, future climatic variables were downscaled under all SSP-RCP scenarios (SSP1-1.9, SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP4-3.4, SSP4-6.0, SSP5-3.4OS, and SSP5-8.5) from MIROC6 and MRI-ESM-2.0, respectively. Downscaled projections based on CMIP6 models are likely to generate the warmer and dryer climate over Ishikari River basin. Increased temperature and less precipitation changes in the far-future period (2070s) may be larger than those in the middle period (2040s). Hence, our results forecast plausible future climate change. The datasets of climatic variables established in this study can be utilized in regional or local hydrologic and environmental modeling, as well as in analyzing the sustainability of ecosystems in further research.